NUMERICAL COMPUTATIONS OF VISCOUS, INCOMPRESSIBLE FLOW 
PROBLEMS USING A TWO-LEVEL FINITE ELEMENT METHOD 
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Abstract. We consider two-level finite element discretization methods for the stream function formulation of the 
Navier-Stokes equations. The two-level method consists of solving a small nonlinear system on the coarse mesh, then 
solving a linear system on the fine mesh. The basic result states that the errors between the coarse and fine meshes 
are related superlinearly. This paper demonstrates that the two-level method can be implemented to approximate 
efficiently solutions to the Navier-Stokes equations. Two fluid flow calculations are considered to test problems which 
have a known solution and the driven cavity problem. Stream function contours are displayed showing the main features 
of the flow. 

Key words. Two-level method, Navier-Stokes equations, flnite element, stream function formulation, Reynolds 
number 

AMS subject classifications. 65N35, 76M30, 76D05 

1. Introduction. The numerical treatment of nonlinear problems that arise in areas such as 
fluid mechanics often requires solving large systems of nonlinear equations. Many methods have been 
proposed that attempt to solve these systems efficiently; one such class of methods are two-level 
methods that reduce the computational time. The computational attractions of the methods are that 
they require the solution of only a small system of nonlinear equations on coarse mesh and one linear 
system of equations on fine mesh. Apparently, the two-level method was proposed first in [17, 16, 15] 
and used for semilinear elliptic problems. The method was implemented for the velocity-pressure 
formulation of the Navier-Stokes equations in [11, 12, 13] and for the stream function formulation of 
the Navier-Stokes equations in [8, 18, 7]. 

The Navier-Stokes equations may be solved using either the primitive variable or stream function 
formulation. Here we use the stream function formulation. The attractions of the stream function 
formulation are that the incompressibility constraint is automatically satisfied, the pressure is not 
present in the weak form, and there is only one scalar unknown to solve for. The standard weak 
formulation of the stream function version first appeared in 1979 in [10]. In this direction, Cayco and 
Nicolaides [5, 4, 3] studied a general analysis of convergence for this standard weak formulation. 

The goal of this paper is to demonstrate that the two-level method can be implemented to ap- 
proximate solutions for incompressible viscous flow problems with high Reynolds number. 

2. Governing Equations. Consider the Navier-Stokes equations describing the flow of an in- 
compressible fluid 

(2.1) 
(2.2) 
(2.3) 

where u = (mi, M2) and p denotes the unknown velocity and pressure field, respectively, in a bounded, simply 
connected polygonal domain C it!^ J is a given body force and Re is the Reynolds number. 
The introduction of a stream function ip{x,y) defined by 

dip dip 
oy ox 

means that the continuity equation (2.2) is satisfied identically. The pressure may then be eliminated 
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from (2.1) to give 

Re"^A^V - '^y^ipx + i^x^i'y = curl/ in (2.4) 

V' = on on, (2.5) 

^ = on on, (2.6) 

an 

where n represents the outward unit normal to Q. In order to write (2.4)-(2.6) in a variational form, 
we define the Sobolev spaces 

H\n) = {v:ve L^{n),Dv € L^{n)}, (2.7) 
H^(n) = {v:vGH\n),v = Oon dfl}, (2.8) 

H^{n) = {v:ve L'^in), Dv e L^{n),D'^v e l^(0)}, (2.9) 
H^{n) = {ve H^{n) ■.v = — = o, on on}, (2.10) 

where L^(0) is the space of square integrable functions on ft and D represents differentiation with 
respect to x or y. For each (p G H^{fl), define curl = ^ ^. The standard weak form of 
equations (2.4)-(2.6) is: 

Find V e H^{Cl) such that, for all(p£H^{n), 



where 



a{i), (j)) = Re-^ / Alp- A(j) dCt, 
Jn 

Jn 

l{(f)) = (/,curl(/)) = / / • curl^ dQ.. 



3. Finite Element Discretization. For the standard finite element discretization of (2.11) we 
choose subspace c i?o(ri). We then 

seek e such that, for all (j)" G X", 

One can prove existence and uniqueness for the solution of the discrete problem (3.1) ; see [5] for 

low Reynolds numbers and sec [8, 7] for high Reynolds numbers. 

Once the finite clement spaces arc prescribed, the discrete problem (3.1) reduces to solving a 
system of nonlinear algebraic equations which has a Jacobian which is large, sparse and bounded. 
Various iterative methods can be used to solve the nonlinear problem (3.1). For example, a standard 
approach is to use Newton's method to linearize (3.1) for a fixed mesh of size h. Since only one fixed 
mesh spacing is used, we will refer to this approach as a one-level approach. In the next section we 
will describe a two- level method for solving the discrete problem (3.1). 

4. Two-level Method. We consider the approximate solution of (2.4) by a two-level, finite- 
element procedure. Let X^,X^ c Hq{Q) denote two conforming finite-element meshes with H ^ h. 
The method we consider computes an approximate solution ^p'^ in the finite-element space X^ by 
solving one linear system for the degrees of freedom on X'^. This particular linear problem requires 
the construction of a finite-element space X^ upon a very coarse mesh of width ^ h\ and then the 
solution of a much smaller system of nonlinear equations for an approximation in X^ . The solution 
procedure is then given in Algorithm 1. 



3 



Step 1. Solve the nonlinear system on coarse mesh for G X^: 

a {ij", (j)") + b {ip", = (/,curl (P"^ , for all (j)" G X" . (4.1) 

Step 2. Solve the linear system on fine mesh for ip^ G X^: 

a (i)^, (jJ") +b{^", ^p'\ = (/, curl cj)''^ , for all c/)^ G X^. (4.2) 



Algorithm 1: The Two-Level Algorithm 

The inclusion X^ C Hq{Q) requires the use of finite-element functions that are continuously 
differentiable over Q. We shall give some examples of finite-element spaces for the stream function 
formulation (sec [6, 3]). Wc will impose boundary conditions by setting all the degrees of freedom 
at the boundary nodes to be zero and the normal derivative equal to zero at all vertices and nodes on 
the boundary. 

Argyis Triangle. The functions are quintic polynomials within each triangle and the 21 degrees of 
freedom are chosen to be the function value, the first and second derivatives at the vertices, and the 
normal derivative at the midsides. 

Clough- Tocher Triangle. Here we subdivide each triangle into three triangles by joining the vertices 

to the ccntroid. In each of the smaller triangles, the fmictions are cubic polynomials. There are then 
30 degrees of freedom needed to determine the three different cubic polynomials associated with the 
three triangles. Eighteen of these are used to ensure that, within the big triangle, the functions are 
continuously differentiable. The remaining 12 degrees of freedom arc chosen to be the function values 
and the first derivatives at the vertices and the normal derivative at the midsides. 

Bogner-Fox-Schmit Rectangle. The functions are bicubic polynomials within each rectangle. The 
degrees of freedom are chosen to be the function value, the first derivatives, and the mixed second 
derivative at the vertices. We set the function and the normal derivative values equal to zero at all 
vertices on the boundary. 

Bicubic Spline Rectangle. The functions are the product of cubic splines. These functions are 
bicubic polynomials within each rectangle, are twice continuously differentiable over and their 
degrees of freedom are the function values at the nodes (plus some additional ones on the boundary) . 

The question which automatically arises is how to choose X^ and X'^ so that we obtain opti- 
mal accuracy. This question was addressed from a theoretical standpoint in [8] and the results are 
summarized below. In [8] , it was proven that the algorithm produces an approximate solution which 
satisfies the error bound 

|V-V''|2<c| ,inf^JV'-'!^''|2 + |ln/i|^^^- IV'-V'^lij- (4.3) 

As an example, consider the case of the Clough- Tocher triangle. For this element we have the following 
inequalities: 

<c/i^-^' (j = 0,l,2), 
\i^-i^"\.<CH^-i (j = 0,l,2). 

Thus, if we seek an approximate solution ^p^ with the same asymptotic accuracy as '0'' in | • I2, the 
above error bound shows that the superlinear scaling, between coarse and fine meshes, 

h = 0(^H''/^\lnH\'/^y (4.4) 
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suffices. Analogous scalings between coarse and fine meshes can be calculated from (4.3) by balancing 
error terms on the right-hand side of (4.3) in the same way. For each of the elements described above, 
we give, in Table 4.1, the scaling between coarse and fine meshes. 



Element 


1 V-V'" I2 


1 li 


Scahng 


Argyris triangle 






h 


\nh 


-1/4= 0(if5/2) 


Clough- Tocher triangle 






h 


\nh 


-1/4= 0(7i-3/2) 


Bogner-Fox-Schmit rectangle 






h 


Inh 


-1/4= 0(if3/2) 


Bicubic spline rectangle 






h 


Inh 


-1/4= 0(ff3/2) 



Table 4.1 
Scaling between coarse and fine meshes 



5. Numerical Examples. In this section we describe some numerical results obtained by imple- 
menting the two-level algorithm for which we have an exact solution and the second is the well-known 
driven cavity problem. We chose the later problem because there are numerous results in the literature 
with which to compare. 

For both examples, the region is the unit square {0<a:<l, < y < 1} and for the finite 
element discretization we use the Bogner-Fox-Schmit elements. In order to compare the efficiency of 
the proposed method all linear and nonlinear systems were solved in the same way. All nonlinear 
problems were solved by Newton's method until the norm of the difference in successive iterates and 
the norm of residual were within a fixed tolerance. In each Newton's iteration, we need to solve a 
linear system. The resulting linear system is non-symmetric whose symmetric part is positive definite. 
Moreover, the resulting matrix is sparse matrix. We choose the Bi-Conjugate Gradient Stabilized 
method (BICGSTAB) which requires two matrix-vector products and four inner products in each 
iteration. BICGSTAB is given and discussed in [2]. When solving the linearized problem with a mesh 
spacing h we need the solution generated on a mesh with spacing H. To do this we interpolate 
the solution onto the grid with spacing h. 

Example 1 

We consider as a test example the 2D Navier-Stokes equations (2.1)- (2.3) on the unit square 
V, = (0, 1)2 where we define the right hand side by / := —Re~^Au+ {u- \/)u + \/p with the following 
prescribed exact solution 

with ^(a;,y)=a;2(a;-l)V(2/-l)^ 
p = x^ - 0.5. 

For this test problem, all requirements of the theory concerning the geometry of the domain and the 
smoothness of the data are satisfied. Moreover, the stream function ^/'(a;,y) satisfies the boundary 
conditions of the stream function equation of the Navier-Stokes equations. 

Our goal in this test is to validate the code and the properties and merits of the two-level method 
as compared with the one-level method. In all numerical calculations in this example we have used the 
Bogner-Fox-Schmit elements with Re = 10 and tol = 10-^. We pick three values of h. They are |, 
and jq. The cpu-time, munber of Newton's iterations, number of Bicgstab iterations, the L^-error 
and H'^-eiroT of the stream function for the one-level method for different values of h are tabulated 
in Table 5.1. Table 5.2 shows cup-time, number of Newton's iterations and the number of Bicgstab's 
iterations for each linear solver for the two- level method. Figure 5.1 shows, for fixed Re = 10 using 
the one- level and two- level method, the streamlines for h= 
Remetrks 

1. From Tables 5.1 and 5.2, the cpu-time for the two-level method is much smaller than the 
corresponding cpu-time for the one-level method. For h = hiji, we save about 57%. For 
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example, in /i = ^ we save about 73%. For example, in /i = ^ we need to solve a nonlinear 
system of 1156 equations which requires solving throe linear systems of equations of order 1156. 
The corresponding two-level method requires solving a nonlinear system of 324 equations and 
a linear system of 1156 equations. We anticipate the savings to increase as the mesh is further 
refined. 

2. From Figures 5.1 , both columns are exactly the same, which means that the two-level method 
produces a solution with the same quality as the one-level method. 

3. From Tables 5.1 and 5.2, both iJ^-error and iJ^-error are of the same order, which means 
that the velocity field is of the same error and quality in both methods since u = tpy and 
V = -V'x- 
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cpu time 
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W^-i'" h,h 
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69.74 
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26,30,24 


3.94x10-'"^ 


4.23xl0--^ 


1.41x10-1 


14 


294.35 


3 


50,56,55 


2.41x10-^ 


1.74x10"^ 


8.87x10-^ 


i 

Ifi 


992.23 
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63,75,67 


3.55x10-'^ 


1.41x10"^ 


8.08x10"^ 



Tabic 5.1: One level method, where ni = number of Newton's 
iteration, nb = number of Bicgstab iterations 
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cpu time 


noc 


nof 


u-r \kh 


II V'-^'^ 111,/. 


II h,h 


i i 


29.91 


10,18,15 


50 


1.40x10"'' 


4.09x10""^ 


1.42x10"^ 




128.28 


22,24,21 


118 


2.31x10"* 


1.73x10"^ 


8.70x10"^ 


8' 16 


267.33 


26,30,24 


169 


1.70x10"* 


1.41x10"^ 


7.94x10"^ 



Table 5.2: Two level method, where noc = number of Bicgstab 
iterations in coarse, nof = number of Bicgstab iterations in fine 



Example 2 

In this example, we consider the problem which was described in Example 1. The exact solution 
is very smooth and does not depend on the Reynolds number. The point of these tests is to increase 
the Reynolds number with h fixed and test the robustness of the method. Our goal in this test is to 
determine the validation of the code and the norm behavior when Re is varying. Schieweck [14] tested 
his code with this problem. He used two nonconforming finite element approximations of upwind type 
for the velocity-pressure formulation. 

The numerical computations of this example were obtained using a Sun Ultra 2 with 2 200 Mhz 
ultrasparc processor running Solaris 2.5.1. In all numerical calculations, we used the Bogncr-Fox- 
Schmit rectangles. The streamlines for Re = 100, 1000 and 2000 were obtained with 17 x 17 grid 
points on the coarse mesh and 33 x 33 grid on the fine mesh. Hence, a mesh of 256 elements and a 
mesh of 1024 elements were used in this test for the case of Bgner-Fox-Schmit rectangles. 

Table 5.3 represents the i^-error, H^-eivoT and iJ^-error. The streamlines are plotted in Figure 5.2 
for Re = 100, 1000, 2000. Figure 5.2 shows that an increase in the Reynolds number will affect the 
streamlines and increase the number of corner contours. 

ExEimple 3 

Cavity flows have been a subject of study for some time. These flows have been widely used as test 
cases for validating incompressible fluid dynamics algorithm. Corner singularities for two-dimensional 
fluid flows are very important since most examples of physical interest have corners. For example, 
singularities of most elliptic problems develop when the boundary contour is not smooth. In this 
example, we consider the driven flow in a rectangular cavity when the top surface moves with a 
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constant velocity along its length. The upper corners where the moving surface meets the stationary 
walls are singular points of the flow at which the horizontal velocity is multi- valued. The lower corners 
are also weakly singular points. 

We consider a domain CI = [0,1] x [0,1] with no-slip boundary conditions, i.e., u = u = in 
all boundaries except y = 1, where u = 1. This problem has been studied and addressed by many 
researchers including Ghia, Ghia, Shin [9], and J.E. Akin [1]. The numerical computation of this 
example was obtained using a Sun Ultra 2 with 2 200 MHz ultrasparc processor running Solaris 2.5.1. 
Bogner-Fox-Schmit elements are used with 17 x 17 grid points on the coarse mesh and 33 x 33 grid 
points on the fine mesh. The streamlines for Re = 1, 10, 50, 100 are plotted in Figure 5.3. Figure 5.4 
shows the w-velocity lines through the vertical line x = 0.5 and t;-velocity lines through the horizontal 
line y = 0.5. 
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II ^-'A'^ 111,/. 


11^-^" II2,/. 
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\ Ifi ' ,32 ) 


61,66,43 
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V 16 ' 32 / 
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4.44x10-"* 


5.00x10"^ 


5.11x10"^ 


200 


^) 

32' 


53,254,421, 
668,1287,454 


4356 


3.06x10"* 


4.52x10"^ 


5.07x10"^ 
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(- -) 

V 16 ' 32 
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9.62x10"* 
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1.93 xlO"^ 
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4.93 xlO"i 



Table 5.3: Two level for test problem, where noc = number of 
Bicgstab iterations in coarse , nof = number of Bicgstab iteration 
in fine 
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Fig. 5.1. Left-side: streamlines for h = with Re = 10 using the one-level method. Right-side: streamlines 

for (H^h) = (^,^) '(7)]^) '(^'1^) '^'^^^ ^6 = 10 using the two-level method. 



Fig. 5.4. Cavity Problem :(above)u-velocity lines through the vertical line x 
and Ghia lines for Re=400, (below) v-velocity lines through the horizontal line y 
and Ghia lines for Re=400. (courtesy U. Ghia [9]) 



= 0.5 for different Reynolds number 
= 0.5 for different Reynolds number 
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